Global changes have affect a tremendous effect communities
Understanding assembly and disassembly rules
Give a more complex and complete picture of what is going on ?
Many papers report temporal trends of local species richness and turnover
filtered_dataset <- map(filtered_dataset,
~.x %>%
filter(siteid %in% unique(modelling_data$siteid))
)
uabun <- tabyl_df(filtered_dataset$site_quali, "unitabundance") %$%
setNames(.$percent, .$unitabundance)
protocol <- tabyl_df(filtered_dataset$site_quali, "protocol")
electro <- str_detect(filtered_dataset$site_quali$protocol, "Electro") %>%
sum / nrow(filtered_dataset$site_quali) * 100
year_samp <- modelling_data %>%
summarise(enframe(summary_distribution(year))) %>%
mutate(value = round(value)) %>%
deframe()
year_samp_hft <- modelling_data %>%
summarise(
be93 = sum(year <= 1993) / n(),
a93 = sum(year > 1993) / n(),
bw9309 = sum(year > 1993 & year <= 2009) / n(),
a09 = sum(year > 2009) / n()
) %>%
mutate(across(everything(), scales::percent)) %>%
pivot_longer(everything(), names_to = "name", values_to = "value") %>%
deframe()
Unique protocol and unique unit of abundance by site
One sampling by year and site
In each site, samplings should be in the median (more or less one month and half, i.e. 1.5 month)
5 samplings minimum by site, span 10 years
5359 sites
98.1% of the sites were monitored with electrofishing
The year of samplings goes from 1955 to 2019 with first, second and third quartile being respectively 1999, 2006, 2012.
loc <- filtered_dataset$location %>%
left_join(filtered_dataset$site_quali, by = "siteid") %>%
st_as_sf(coords = c("longitude", "latitude"),
crs = 4326)
fig_sites <- paste0(
"Distribution of the sites.")
ggplot(data = world) +
geom_sf() +
geom_sf(data = loc, aes(color = protocol), shape = 1) +
# coord_sf(xlim = bb[c("xmin", "xmax")],
# ylim = bb[c("ymin", "ymax")],
# datum = 4326, expand = TRUE) +
theme(legend.position = "bottom") +
labs(title = paste0("Number of sites: ", nrow(loc)))
Figure 3.1: Distribution of the sites.
Total abundance was reported in four different units: Count, Count by 100m^2, Leslie index & CPUE. The two first units account for 80% of the data.
We computed species richness with coverage-based correction proposed by (???). The Chao index correct the species richness for the sampling coverage. The sampling coverage of each sampling is evaluated by the total number of individual and the number of singleton and (double) species. All the sampling events were set at the same coverage of .80, interpolation and rarefaction for respectively low and high coverage sampling. Chao index was logged to express the temporal trends of species richness of percentage variation. The temporal trends estimated from uncorrected species richness gave the same results (Figure SXX).
Betadiversity metrics compares the composition of two communities. In our case, we compare the composition of each time step relatively to the baseline, i.e. the first year of sampling.
Turnover metrics based on presence/absence can be decomposed in appearance/disappearance and nestedness / turnover.
appearance: \(\dfrac{S_{imm}}{S_{tot}}\), disappearance: \(\dfrac{S_{lost}}{S_{tot}}\)
Turnover: \(J_t = \dfrac{2 * min(S_{lost}, S_{imm})}{S_{common} + (2 * min(S_{lost}, S_{imm})}\)
Nestedness: \(SER_r - J_t\)
var_pca <- dimnames(pca_riv_str$rotated$loadings)[[1]]
clean_var_pca <- get_river_atlas_significant_var()[
get_river_atlas_significant_var() %in% var_pca]
high_load_rc1 <- pca_riv_str$rotated$loadings[
pca_riv_str$rotated$loadings[, "RC1"] > .90,
"RC1"]
clean_high_load_rc1_var <- get_river_atlas_significant_var()[
get_river_atlas_significant_var() %in% names(high_load_rc1)]
fig_cap_pca <- "Rotated PCA on stream structure. The first axis is related to
the distance from the source, stralher order and discharge (m3/s). The second
axis is related to the altitude and slope degree of the site. Only the first
axis was included in the analysis"
my_pca_plot(pca_riv_str$rotated)
Figure 4.1: Rotated PCA on stream structure. The first axis is related to the distance from the source, stralher order and discharge (m3/s). The second axis is related to the altitude and slope degree of the site. Only the first axis was included in the analysis
h_hft_cap <- "Distribution of log base 2 ratio between human footprint in 1993 and
2009. -1 and +1 means an Human footprint respectively divided and multiplied by
two"
modelling_data %>%
group_by(siteid) %>%
summarise(hft_ix_c9309_log2_ratio = unique(hft_ix_c9309_log2_ratio)) %>%
ggplot(aes(x = hft_ix_c9309_log2_ratio)) +
geom_histogram() +
labs(
x = "Human footprint ratio (1993-2009, Log base 2)",
y = "Number of sites"
) +
theme_cowplot()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 4.2: Distribution of log base 2 ratio between human footprint in 1993 and 2009. -1 and +1 means an Human footprint respectively divided and multiplied by two
s_hft_diff <- round(summary_distribution(modelling_data$hft_ix_c9309_log2_ratio, na.rm = T), 2)
In each site and at each sampling point, we computed the dissimilarity (Jaccard, Appearance, Disappearance, Turnover & Nestedness) between each time step and the reference time, i.e. the first year of sampling. Regarding the reference time step, there are two different strategies. First is to include it, i.e. setting a dissimilarity of 0 at each beginning of a time series (Blowes et al. 2019). Another strategy is to exclude the first year from the analysis (Dornelas et al. 2014). The first option logically constrains the analysis toward increasing dissimilarity as the second can capture recovery from perturbation soon after the first samplings. I guess that the first option makes more sense more it removes to sensibility to community changes happening just after the references year.
However as the relationship between most dissimilarity indices and time is highly non linear (like a Michaelis-Menten), we logged the number of year since the first year of sampling (plus one) to linearize the relationship. AIC of the models with logged years were up to twice lower (Table SXX).
Jaccard dissimilarity, nestedness, turnover, appearance, disappearance and hillebrand were modelled in two ways, Beta and gaussian regression. Beta regressions are more appropriate for distribution bounded by [0,1] and are linearized by a logit link. Logit links have the cons that coefficients are less easy to interpret. So we run both gaussian and beta regressions and investigate the sensitivity of our results to modelling choice. Values of turnover were slightly transformed such as values of one and zero fitted beta regression requirements (\(x' = \dfrac{x \times (N - 1) + 0.5}{N}\)).
The model is defined like this:
\[ \begin{align} \\ B_i &= \beta_0Y_i + \beta_1Y_iR_i + \beta_2Y_iH_i + \beta_3Y_i{HD}_i + \\& \beta_4Y_iR_iH_i + \beta_5Y_iR_i{HD}_i + \\& \beta_6Y_iH_i{HD}_i + \epsilon \end{align} \]
\[ \begin{align} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0i} \end{align} \]
"biodiv_facet ~
0 + year/river_structure + year/human_footprint_diff +
(0 + year | basin/site)"
#> [1] "biodiv_facet ~ \n0 + year/river_structure + year/human_footprint_diff +\n(0 + year | basin/site)"
year: logged number of year (+1) since the beginning of the monitoring (site scale)river_structure: PCA from upstream to downstream (centered and scaled)human_footprint_diff: difference between 1993 and 2009 (centered and scaled)basin: Hydrographic basinsite: siteAs temporal trends of biodiversity facets had a Michaelis-Menten shape, we linearized the trends by taking the log of year number added to 1. Models with logged year had a lower AIC, up to twice (Table SXX).
At each site, all the dissimilarity indices had a value of 0 the first year. Then we constrained the intercept to 0 and kept only random effect on the slope. Models comparison shows that the results are qualitatively the same (Appendix SXX).
At the difference with dissimilarity index, we modelled species richness with the intercept and the main effect of ecological drivers. We modelled logged species richness with Chao index.
\[ \begin{align} \\ B_i &= \alpha + \beta_0Y_i + \beta_1R_i + \beta_2Y_iR_i + \beta_3H_i + \beta_4Y_iH_i + \\& \beta_5Y_iR_iH_i + \beta_6Y_iR_i{HD}_i + \\& \beta_7Y_iH_i{HD}_i + \epsilon \end{align} \]
\[ \begin{align} \\ \alpha&=\alpha_0 + u_{0i|j} + u_{0j} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0j} \end{align} \]
In the dataset, total abundance was described as raw count, counts by \(100m^2\), CPUE, and Leslie index representing respectively XX% of the sites. To account for this hererogeneity, we included the unit of measurement as a main affect explaining total abundance but also as an interaction with the year effect, in a similar way than previous studies (Klink et al. 2020). Similar to species richness modelling, we modelled the log of total abundances.
\[ \begin{align} \\ B_i &= \alpha + \beta_0Y_i + \beta_1U_i + \beta_2Y_iU_i + \beta_3R_i + \\& \beta_4Y_iR_i + \beta_5H_i + \beta_6Y_iH_i + \\& \beta_7Y_iR_iH_i + \beta_8Y_iR_i{HD}_i + \\& \beta_9Y_iH_i{HD}_i + \epsilon \end{align} \]
where \[ \begin{align} \\ \alpha&=\alpha_0 + u_{0i|j} + u_{0j} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0j} \end{align} \]
All the statistical models were computed in R with the glmmTMB package. In
order to facilitate the comparison of the effects, the reported slope
coefficients were standardised as \(r_\delta = \beta \dfrac{\s_x}{\s_y}\), \(\beta\)
being the unstandardised slope, \(s_x\) and \(s_y\) respectively the standard
deviation of the explicative and response variable.
The fixed effect coefficients for each level of random effects were computed as
sum of the random effects and its corresponding fixed effects, known as BLUPs
(Best Linear Unbiased Predictions, coef.merMod in R). However, the associated
standard errors and confidence intervals are not computable on BLUPs. As an
example, it means that while we are able to compute the temporal trends as each
site and basin, those temporal trends have no associated errors.
tar_load(site_no_drivers_inla_tot)
The estimated temporal trends of total abundance were corrected according to abundance units, with the estimated fixed effect from the model.
kmeans_opt_nb_gr_sil_cap <- "Objective function value according to the
proportion of outliers data trimmed and the number of cluster(colored lines). Up
to 6 cluster, the proportion of trimmed data does not change the ranking of the clustering. We chose to stick with 6 clusters with a removal 0.05 of the outliers"
kmeans_sites <- "Distribution of the sites in discriminant space, colored by
their cluster assignement (left). Stilhouette profile. Black points represent sites that have been trimmed
at the 5% threshold."
Site were clustered based on the temporal trends of biodiversity facets using
trimmed robust cluster method (???). Temporal trends were scaled but not
centered before clustering. According to the optimal function, the optimal
number of clusters was between 6 and 8 (Fig. ??). We used 25 starting
positions and a maximum of 100 iterations. We trimmed 5% of outliers in the
multidimentional space. We set cluster assignement as NA for sites that had an
assignement improvement of 50% in another cluster. Clustering was computed with
tclust R package.
restricted_dimension <- c(
"hillebrand_dis_scaled", "log_total_abundance",
"log_chao_richness", "perc_exo_sp", "perc_exo_abun")
tu <- rbind( gaussian_inla_std_effects,
gaussian_inla_exo_std_effects) %>%
filter(!str_detect(term, "unitabundance|Intercept")) %>%
filter(response %in% restricted_dimension) %>%
mutate(
term = str_replace_all(term, get_model_term_replacement()),
response = get_var_replacement()[response]
)
library(ggbreak)
p1 <- tu %>%
filter(facet == "main") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect() +
scale_x_break(c(.25, .65))
p2 <- tu %>%
filter(facet == "interaction") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect() +
scale_x_break(c(.10, .17))
p3 <- tu %>%
filter(facet == "dbl_interaction") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect(.,
legend_present = TRUE,
xaxis_title = TRUE
) +
labs(x = "Standardized coefficients")
leg_dis <- get_legend(p3 +
guides(colour = guide_legend(nrow = 2) ) +
theme(legend.position = "bottom")
)
fig1 <- plot_grid(
print(p1),
print(p2),
p3 + theme(legend.position = "none"),
leg_dis,
rel_heights = c(1, 1, 1, .3),
ncol = 1)
fig1
tar_load(site_no_drivers_inla)
cor_site_tps <- cor(site_no_drivers_inla,
method = "spearman"
)
dimnames(cor_site_tps) <- map(
dimnames(cor_site_tps),
~get_var_replacement()[.x]
)
#corrplot::corrplot(cor_clust, type = "upper", diag = FALSE, tl.srt = 15)
cor_p <- ggcorrplot::ggcorrplot(
cor_site_tps, hc.order = TRUE, type = "lower",
lab = TRUE, tl.srt = 15, legend.title = "Spearman Corr") +
theme(legend.position = c(.45, .95), legend.direction = "horizontal") +
guides(colour = guide_legend(title.position = "top"))
mod_slp_std_cap <- paste0("(A) Estimated fixed effects of the gaussian models.
Coefficients were standardized by refiting the models with scaled response and
explicative variables. Bars represent CI at 95% computed with Wald
approximation. (B) Interactions effects of pressures on temporal trends
(Antagonistic or synergetic).")
color_cl_scale <- setNames(pal2, unique(site_cl_rm$cl))
pca_clust_cap <- "PCA biplot on site temporal trends (model without
drivers). The sites are colored by their cluster belongings."
tar_load(c(pca_clust, site_cl_rm))
pca_clust <- compute_rotated_pca(site_no_drivers_inla_tot, naxis = 5)
p_rotated_pca_clust <- my_pca_plot(
pca_clust$rotated,
replace_var = get_var_replacement())
tar_load(p_pca_clust)
p_pca_clust <- plot_pca_clust(
.data = pca_clust$rotated,
site_cl = site_cl_rm,
replace_var = get_var_replacement(),
size_arrows_segment = 1,
label_size = 2.5,
alpha_point = .2,
lim_x_y = c(-3.5, 3.5),
force_pull = 1,
force = 80
)
axis_l <- list(
c("RC1", "RC2"),
c("RC1", "RC3"),
c("RC1", "RC4"),
c("RC1", "RC5"),
c("RC2", "RC3"),
c("RC2", "RC4"),
c("RC2", "RC5"),
c("RC3", "RC4"),
c("RC3", "RC5"),
c("RC4", "RC5")
)
p_pca_cl_ell <- map(axis_l,
~plot_pca_clust(
.data = pca_clust$rotated,
site_cl = site_cl_rm,
add_point = TRUE,
add_ellipse = TRUE,
xaxis = .x[1], yaxis = .x[2],
ctb_thld = .2,
replace_var = get_var_replacement(),
size_arrows_segment = 1,
label_size = 2.5,
alpha_point = .2,
lim_x_y = c(-3.5, 3.5),
force_pull = 1,
force = 80,
var_scaling_factor = 3.5
) +
theme(legend.position = "bottom") +
scale_color_manual(values = color_cl_scale)
)
leg_cl <- get_legend(p_pca_cl_ell[[1]])
p_pca_cl_ell <- map(p_pca_cl_ell, ~.x + theme(legend.position = "none"))
p_pca_cl_top <- plot_grid(plotlist = p_pca_cl_ell[c(1, 5, 6, 7)], nrow = 2)
#> Warning: Removed 55 rows containing non-finite values (stat_ellipse).
#> Warning: Removed 55 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_ellipse).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_ellipse).
#> Warning: Removed 1 rows containing missing values (geom_point).
p_pca_cl_bt <- target_bp_cl_dist(cl_obj = site_cl_rm) +
scale_fill_manual(values = color_cl_scale) +
scale_color_manual(values = color_cl_scale) +
theme(axis.text.x = element_text(angle = 20, vjust = 0))
plot_grid(
p_pca_cl_top,
p_pca_cl_bt,
nrow = 2,
rel_heights = c(1, .5)
)
Figure 6.1: PCA biplot on site temporal trends (model without drivers). The sites are colored by their cluster belongings.
site_cl_rm_loc <- site_cl_rm %>%
select(siteid, cl) %>%
left_join(loc, by = "siteid") %>%
filter(ecoregion %in% c("Palearctic", "Nearctic", "Australasia")) %>%
left_join(
measurement_exo %>%
distinct(siteid, basin_name),
by = "siteid")
n_clust_basin <- site_cl_rm_loc %>%
group_by(ecoregion, country, basin_name) %>%
summarise(
n = n(),
div = length(unique(cl)),
evenness = ifelse(div == 1, 1, diversity(table(cl)) / log(div))
)
#> `summarise()` has grouped output by 'ecoregion', 'country'. You can
#> override using the `.groups` argument.
basin_name_palearctic <- n_clust_basin %>%
filter(ecoregion == "Palearctic", div == 6) %>%
arrange(desc(n), desc(evenness))
#Derwent.UK, Loire, Seine
fq_pal_bas <- site_cl_rm_loc %>%
filter(basin_name %in% c("Derwent.UK", "Umealven", "Seine")) %>%
group_by(basin_name, cl) %>%
summarise(freq = n()) %>%
mutate(pc = prop.table(freq)) %>%
ungroup()
#> `summarise()` has grouped output by 'basin_name'. You can override using
#> the `.groups` argument.
fq_pal <- site_cl_rm_loc %>%
filter(ecoregion == "Palearctic") %>%
group_by(cl) %>%
summarise(freq = n()) %>%
mutate(pc = prop.table(freq), basin_name = toupper("Palearctic")) %>%
ungroup()
fq_tab <- rbind(
fq_pal, fq_pal_bas
)
get_fq_tab_from_basin <- function(
x = site_cl_rm_loc,
ecoregion = NULL,
basins = NULL) {
fq_pal_bas <- x %>%
filter(basin_name %in% basins) %>%
group_by(basin_name, cl) %>%
summarise(freq = n()) %>%
mutate(pc = prop.table(freq)) %>%
ungroup()
ecoregion2 <- ecoregion# to not mess up filter
fq_pal <- site_cl_rm_loc %>%
filter(ecoregion == ecoregion2) %>%
group_by(cl) %>%
summarise(freq = n()) %>%
mutate(
pc = prop.table(freq),
basin_name = toupper(ecoregion)
) %>%
ungroup()
fq_tab <- rbind(fq_pal, fq_pal_bas) %>%
mutate(
basin_name = factor(
basin_name, levels = c(toupper(ecoregion), basins)),
ecoregion = ecoregion
)
return(fq_tab)
}
basin_name_palearctic <- n_clust_basin %>%
filter(ecoregion == "Nearctic", div >= 5) %>%
arrange(desc(n), desc(evenness))
tab_palearctic <- get_fq_tab_from_basin(
x = site_cl_rm_loc,
ecoregion = "Palearctic",
basins = c("Seine", "Derwent.UK", "Umealven")
)
#> `summarise()` has grouped output by 'basin_name'. You can override using
#> the `.groups` argument.
tab_nearctic <- get_fq_tab_from_basin(
x = site_cl_rm_loc,
ecoregion = "Nearctic",
basins = c("Mississippi", "Saint.Laurent", "Columbia")
)
#> `summarise()` has grouped output by 'basin_name'. You can override using
#> the `.groups` argument.
basin_name_australasia <- n_clust_basin %>%
filter(ecoregion == "Australasia", div >= 5) %>%
arrange(desc(n), desc(evenness))
tab_australasia <- get_fq_tab_from_basin(
x = site_cl_rm_loc,
ecoregion = "Australasia",
basins = c("Murray.Darling", "Brisbane", "Logan")
)
#> `summarise()` has grouped output by 'basin_name'. You can override using
#> the `.groups` argument.
tot_tab <- rbind(
tab_palearctic,
tab_nearctic,
tab_australasia) %>%
group_by(basin_name) %>%
mutate(
total_n = sum(freq)) %>%
ungroup() %>%
arrange(ecoregion, total_n) %>%
mutate(
basin_name_n = paste0(basin_name, " (", total_n, ")"),
basin_name_f = factor(basin_name_n, levels = unique(basin_name_n))
)
p_tot_tab <- tot_tab %>%
ggplot(aes(x = pc, y = basin_name_f, fill = as.factor(cl))) +
geom_col() +
labs(x = "Proportion") +
scale_x_continuous(expand = c(0, 0)) +
theme_cowplot() +
theme(axis.title.y = element_blank(), legend.position = "none") +
facet_grid(rows = vars(ecoregion), scale = "free_y") +
scale_fill_manual(values = color_cl_scale)
tar_load(rivers10)
world_site_sf$cl <- world_site_sf$site %>%
left_join(site_cl_rm[, c("siteid", "cl")], by = "siteid") %>%
mutate(cl = relevel(as.factor(cl), ref = 1))
main_map <- ggplot(world_site_sf$world) +
geom_sf(color = NA) +
theme_void() +
geom_sf(data = rivers10, colour = "lightblue", alpha = .4) +
geom_sf(data = world_site_sf$cl %>%
filter(!is.na(cl)),
aes(fill = cl), pch = 21, color = "white", alpha = .3) +
scale_fill_manual(values = color_cl_scale) +
theme(legend.position = "none")
main_map4zoom <- ggplot(world_site_sf$world) +
geom_sf(color = NA) +
theme_void() +
geom_sf(data = rivers10, colour = "lightblue", alpha = .4) +
geom_sf(data = world_site_sf$cl %>%
filter(!is.na(cl)),
aes(fill = cl), pch = 21, color = "white", alpha = .3, size = 5) +
scale_fill_manual(values = color_cl_scale) +
theme(legend.position = "none")
zoom_map <- list(
north_am = "USA",
eng = "GBR",
south_eu = c("FRA", "ESP", "BEL"),
north_eu = c("FIN", "NOR", "SWE"),
aus = "AUS"
)
zoom_map_coord <- map(zoom_map,
~world_site_sf$cl %>%
filter(country %in% .x) %>%
st_bbox()
)
zoom_map_coord_df <- map_dfr(zoom_map_coord,
~setNames(as.numeric(.x), names(zoom_map_coord[[1]])) %>%
enframe(.),
.id = "zoom") %>%
pivot_wider(names_from = "name", values_from = "value") %>%
mutate(
label = LETTERS[seq_len(length(zoom_map_coord))],
)
main_map_rect <- main_map +
geom_text(
data = zoom_map_coord_df,
aes(x = xmax, y = ymax, label = label),
vjust = 0.5, hjust = 0, size = 7) +
coord_sf(expand = FALSE) +
theme(legend.position = "none")
save_plot(
here("doc", "fig", "global_map.png"),
main_map_rect,
base_height = 3.71,
base_asp = 2.618,
)
p_zoom_map <- map2(zoom_map_coord, zoom_map_coord_df$label,
~main_map4zoom +
coord_sf(
xlim = .x[c("xmin", "xmax")],
ylim = .x[c("ymin", "ymax")],
expand = FALSE) +
annotate("text", x = .x["xmin"], y = .x["ymax"],
label = .y, size = 15, vjust = 1, hjust = 0) +
theme(legend.position = "none")
)
map_cl_cap <- "Distribution of the clusters of biodiversity trends (without
cluster being NAs). Bottom boxplots display the distribution of the biodiversity
trends (scaled) by cluster and facets. On the left is displayed the proportion
of cluster by realm for the three main ones.
"
fig3_left <- p_tot_tab
fig3_right <- plot_grid(
main_map_rect, p_zoom_map$eng, p_zoom_map$north_eu,
p_zoom_map$north_am, p_zoom_map$south_eu, p_zoom_map$aus,
ncol = 3, axis = "l"
)
fig3 <- plot_grid(fig3_left, fig3_right, ncol = 2, rel_widths = c(.5, 1))
fig3
Figure 6.2: Distribution of the clusters of biodiversity trends (without cluster being NAs). Bottom boxplots display the distribution of the biodiversity trends (scaled) by cluster and facets. On the left is displayed the proportion of cluster by realm for the three main ones.
save_plot(
here("doc", "fig", "fig3_left.png"),
fig3_left,
base_height = 3.71,
base_asp = 2.618,
)
save_plot(
here("doc", "fig", "fig3_right.png"),
fig3_right,
base_height = 3.71,
base_asp = 2.618,
)
save_plot(here("doc/fig/fig3.png"),
fig3, base_height = 4, nrow = 1, ncol = 2)
knitr::include_graphics(here("doc/fig/map_cluster.png"))
Blowes, Shane A., Sarah R. Supp, Laura H. Antão, Amanda Bates, Helge Bruelheide, Jonathan M. Chase, Faye Moyes, et al. 2019. “The Geography of Biodiversity Change in Marine and Terrestrial Assemblages.” Science 366 (6463): 339–45. https://doi.org/10.1126/science.aaw1620.
Dornelas, Maria, Nicholas J. Gotelli, Brian McGill, Hideyasu Shimadzu, Faye Moyes, Caya Sievers, and Anne E. Magurran. 2014. “Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss.” Science 344 (6181): 296–99. https://doi.org/10.1126/science.1248484.
Hillebrand, Helmut, Bernd Blasius, Elizabeth T. Borer, Jonathan M. Chase, John A. Downing, Britas Klemens Eriksson, Christopher T. Filstrup, et al. 2018. “Biodiversity Change Is Uncoupled from Species Richness Trends: Consequences for Conservation and Monitoring.” Journal of Applied Ecology 55 (1): 169–84. https://doi.org/https://doi.org/10.1111/1365-2664.12959.
Klink, Roel van, Diana E. Bowler, Konstantin B. Gongalsky, Ann B. Swengel, Alessandro Gentile, and Jonathan M. Chase. 2020. “Meta-Analysis Reveals Declines in Terrestrial but Increases in Freshwater Insect Abundances.” Science 368 (6489): 417–20. https://doi.org/10.1126/science.aax9931.
Vellend, Mark, Lander Baeten, Isla H. Myers-Smith, Sarah C. Elmendorf, Robin Beauséjour, Carissa D. Brown, Pieter De Frenne, Kris Verheyen, and Sonja Wipf. 2013. “Global Meta-Analysis Reveals No Net Change in Local-Scale Plant Biodiversity over Time.” Proceedings of the National Academy of Sciences 110 (48): 19456–9.
tu_full <- rbind(gaussian_inla_std_effects,
gaussian_inla_exo_std_effects) %>%
filter(!str_detect(term, "unitabundance|Intercept")) %>%
mutate(
term = str_replace_all(term, get_model_term_replacement()),
response = get_var_replacement()[response]
)
p1_f <- tu_full %>%
filter(facet == "main") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect() +
scale_x_break(c(.21, .5))
p2_f <- tu_full %>%
filter(facet == "interaction") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect()
p3_f <- tu_full %>%
filter(facet == "dbl_interaction") %>%
arrange(desc(width_bar)) %>%
plot_inla_fixed_effect(.,
legend_present = TRUE,
xaxis_title = TRUE
) +
labs(x = "Standardized coefficients")
leg_dis_f <- get_legend(p3_f +
theme(legend.position = "bottom") +
guides(color = guide_legend(nrow = 4))
)
figS1 <- plot_grid(
print(p1_f),
print(p2_f),
p3_f + theme(legend.position = "none"),
leg_dis_f,
rel_heights = c(1, 1, 1, .2),
ncol = 1)
figS1
Reproducibility receipt
## datetime
Sys.time()
#> [1] "2022-04-28 11:49:01 CDT"
## repository
if(requireNamespace('git2r', quietly = TRUE)) {
git2r::repository()
} else {
c(
system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
system2("git", args = c("remote", "-v"), stdout = TRUE)
)
}
#> [1] "commit b403252aed08d455fd9f7f520a45c81f3fbf0fa4"
#> [2] "Merge: 9018f4f 689d6f1"
#> [3] "Author: Danet <ahdanet@ilstu.edu>"
#> [4] "Date: Wed Apr 27 23:45:10 2022 -0500"
#> [5] ""
#> [6] " merge"
#> [7] " "
#> [8] " Merge branch 'main' of https://github.com/alaindanet/RivFishTimeBiodiversityFacets"
#> [9] " "
#> [10] " # Conflicts:"
#> [11] " # doc/xx-meeting-report.html"
#> [12] "origin\tgit@github.com:alaindanet/RivFishTimeBiodiversityFacets.git (fetch)"
#> [13] "origin\tgit@github.com:alaindanet/RivFishTimeBiodiversityFacets.git (push)"
## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 11 (bullseye)
#>
#> Matrix products: default
#> BLAS: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#>
#> locale:
#> [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
#> [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
#> [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] parallel stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggbreak_0.0.9 hrbrthemes_0.8.0
#> [3] factoextra_1.0.7 ade4_1.7-16
#> [5] tclust_1.4-2 ggeffects_1.0.1
#> [7] report_0.5.0.9000 see_0.6.8.1
#> [9] correlation_0.8.0 modelbased_0.7.1.1
#> [11] effectsize_0.6.0.2 parameters_0.16.0
#> [13] performance_0.8.0 bayestestR_0.11.5
#> [15] datawizard_0.2.3 insight_0.15.0
#> [17] easystats_0.4.3 glmmTMB_1.1.2.9000
#> [19] inlatools_0.0.1.9001 INLA_21.11.22
#> [21] sp_1.4-6 foreach_1.5.1
#> [23] Matrix_1.3-2 future_1.24.0
#> [25] slider_0.2.2 vegan_2.5-7
#> [27] lattice_0.20-41 permute_0.9-5
#> [29] codyn_2.0.5 janitor_2.1.0
#> [31] viridis_0.6.2 viridisLite_0.4.0
#> [33] cowplot_1.1.1 terra_1.5-21
#> [35] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0
#> [37] sf_1.0-6 scales_1.1.1
#> [39] kableExtra_1.3.4.9000 here_1.0.1
#> [41] lubridate_1.8.0 magrittr_2.0.2
#> [43] forcats_0.5.1 stringr_1.4.0
#> [45] dplyr_1.0.8 purrr_0.3.4
#> [47] readr_2.1.1 tidyr_1.2.0
#> [49] tibble_3.1.6 ggplot2_3.3.5
#> [51] tidyverse_1.3.1 tarchetypes_0.3.2
#> [53] targets_0.8.1 conflicted_1.1.0
#> [55] rmarkdown_2.11 nvimcom_0.9-124
#>
#> loaded via a namespace (and not attached):
#> [1] utf8_1.2.2 tidyselect_1.1.1 lme4_1.1-26
#> [4] grid_4.0.5 munsell_0.5.0 codetools_0.2-18
#> [7] units_0.8-0 statmod_1.4.35 withr_2.4.3
#> [10] colorspace_2.0-0 highr_0.9 knitr_1.38
#> [13] rstudioapi_0.13 Rttf2pt1_1.3.8 listenv_0.8.0
#> [16] labeling_0.4.2 emmeans_1.7.2 mnormt_2.0.2
#> [19] polyclip_1.10-0 farver_2.0.3 rprojroot_2.0.2
#> [22] coda_0.19-4 parallelly_1.30.0 vctrs_0.3.8
#> [25] generics_0.1.0 TH.data_1.0-10 xfun_0.30
#> [28] R6_2.5.1 gridGraphics_0.5-1 cachem_1.0.4
#> [31] assertthat_0.2.1 multcomp_1.4-16 gtable_0.3.0
#> [34] globals_0.14.0 processx_3.5.2 sandwich_3.0-0
#> [37] rlang_1.0.2 systemfonts_1.0.0 splines_4.0.5
#> [40] extrafontdb_1.0 TMB_1.7.20 broom_0.7.12
#> [43] reshape2_1.4.4 yaml_2.2.1 modelr_0.1.8
#> [46] backports_1.2.1 extrafont_0.17 tools_4.0.5
#> [49] psych_2.0.12 bookdown_0.24 ggplotify_0.1.0
#> [52] ellipsis_0.3.2 jquerylib_0.1.3 plyr_1.8.6
#> [55] Rcpp_1.0.6 progress_1.2.2 classInt_0.4-3
#> [58] ps_1.6.0 prettyunits_1.1.1 zoo_1.8-8
#> [61] haven_2.5.0 ggrepel_0.9.1 cluster_2.1.1
#> [64] fs_1.5.2 data.table_1.13.6 warp_0.2.0
#> [67] reprex_2.0.1 tmvnsim_1.0-2 mvtnorm_1.1-1
#> [70] patchwork_1.1.1 hms_1.1.1 evaluate_0.15
#> [73] xtable_1.8-4 readxl_1.3.1 gridExtra_2.3
#> [76] compiler_4.0.5 KernSmooth_2.23-18 crayon_1.4.2
#> [79] minqa_1.2.4 htmltools_0.5.2 ggfun_0.0.5
#> [82] mgcv_1.8-34 tzdb_0.2.0 aplot_0.1.2
#> [85] DBI_1.1.1 tweenr_1.0.1 sjlabelled_1.1.7
#> [88] dbplyr_2.1.1 MASS_7.3-53.1 boot_1.3-27
#> [91] cli_3.2.0 igraph_1.3.0 pkgconfig_2.0.3
#> [94] numDeriv_2016.8-1.1 xml2_1.3.2 svglite_2.0.0
#> [97] ggcorrplot_0.1.3 bslib_0.2.4 webshot_0.5.2
#> [100] estimability_1.3 rvest_1.0.2 yulab.utils_0.0.4
#> [103] snakecase_0.11.0 callr_3.7.0 digest_0.6.27
#> [106] cellranger_1.1.0 gdtools_0.2.3 nloptr_1.2.2.2
#> [109] lifecycle_1.0.1 nlme_3.1-152 jsonlite_1.7.2
#> [112] fansi_0.5.0 pillar_1.6.4 fastmap_1.1.0
#> [115] httr_1.4.2 survival_3.2-10 glue_1.6.2
#> [118] iterators_1.0.13 ggforce_0.3.2 class_7.3-18
#> [121] stringi_1.7.6 sass_0.4.0 memoise_2.0.0
#> [124] e1071_1.7-4